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ABSTRACT 


A new Petrov-Galerkin finite element formulation has been proposed for transient 
convection-diffusion problems. Most Petrov-Galerkin formulations take into account the 
spatial discretization, and the weighting functions so developed give satisfactory 
solutions for steady state problems. Though these schemes can be used for transient 
problems, there is scope for improvement. The schemes proposed here, which take into 
account temporal as well as spatial discretization, provide improved solutions. 

Electrophoresis, which involves the motion of charged entities under the influence of 
an applied electric field, is governed by equations similiar to those encountered in fluid 
flow problems, i.e., transient convection-diffusion equations. Test problems are solved 
in electrophoresis and fluid flow. The results obtained are satisfactory. It is also expected 
that these schemes, suitably adapted, will improve the numerical solutions of the 
compressible Euler and the Navier-Stokes equations. 
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CHAPTER 1 
INTRODUCTION 


In recent years finite-element methods have been applied to problems of fluid dynamics, 
and various other transport problems. Some of their advantages over other methods are 
the following: ease in modelling intricate geometries and consistent treatment of the 
boundary conditions. However, problems with convection dominated flows have not 
been solved as successfully as problems in solid mechanics, structures, and heat 
conduction. The reason being that symmetric positive-definite operators - which are 
dominant in differential equations governing solid mechanics problems - lead to a "best 
approximation" of the exact solution by the finite-element method. On the other hand, 
problems in fluid mechanics( and some other transport problems ) are governed by both, 
symmetric positive-definite operators and nonsymmetric operators( convection terms). 
Thus, in convection dominated flows, the "best approximation" property of the Galerkin 
finite element is lost and spurious oscillations are observed. 

It is well known that "upwind " differencing in finite-difference schemes may produce 
a wiggle-free solution, but in most cases the side effect is loss of accuracy due to overly 
diffuse solutions. Efforts to improve the results by adding artificial diffusion terms to the 
differential equations may again result in unacceptable levels of amplitude error or 
smearing of discontinuities. This artificial diffusion concept has been the subject of 
substantial controversy and criticism [1,2,3]. Various "upwind" finite-element schemes 
have been proposed, such as those by Christie et al.[4], Heinrich et al. [5], Hughes[6], 
and Hughes and Atkinson[7]. These methods are based on modified weighting 
functions, modified quadrature rules, and variational principles, which yield satisfactory 
results for one-dimensional problems. However in multi-dimensional cases the solutions 
are overly diffuse especially in the presence of transient or source terms. 

It has been demonstrated that the above mentioned problems can be avoided by using 
proper Petrov-Galerkin (PG) schemes [8,9]. As opposed to regular Galerkin 
formulations, PG formulations involve trial and weighting functions from different 
spaces. Such schemes have been applied to various fluid dynamics and convection- 
diffusion problems by Brooks and Hughes [10], Hughes and Tezduyar[ll], Hughes et 
al. [12], and Hughes, Mallet and Franca[13]. These schemes have been developed to 
obtain steady state solutions, and though they have been tested for transient cases, it was 
felt that there was considerable scope for improvement. 

Chapter 2 highlights the development and testing of new weighting functions which 
improve the performance of transient schemes. The underlying philosopy in 
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conceptualizing such functions is the belief that in addition to depending on spatial 
discretization, the weighting functions must also depend on the temporal discretization. 
Johnson et al. [14] have also approached the problem with a similiar objective in mind. 


1.1 APPLICATIONS 


FLOW FIELD SIMULATION Simulation of flow fields has become exceedingly 
important today with the advent of the space age and the increasing demand for 
understanding of other related properties. Although the governing equations ( i.e., the 
full Navier-Stokes equations ) would require tremendous computational resources, 
depending on specific problems, simplifications can be made leading to feasible models. 
In this thesis a two-dimensional, incompressible flow model has been assumed which 
has a wide range of applications. 

The governing equations are the two momentum equations and the continuity 
equation[15] which can be solved numerically [16]. However by further simplifications, 
the number of transient equations can be reduced to one and the number of unknowns to 
two. This is possible by introducing new variables, the vorticity and the stream 
function[17]. This approach, known as the vorticity- stream function approach, has been 
adopted in Chapter 3. 

ELECTROPHORESIS SEPARATION SIMULATION Electrophoresis involves 
motion of dissolved or suspended material under the influence of an applied electric 
field. The entities may be simple ions, complex macromolecules, colloids, or living 
cells. The rate of migration depends on factors such as the amount of charge, the size 
and shape of the entity, and the solvent properties.These properties are conveniently 
combined into a single factor, termed the "particle mobility". Thus, when particles of 
different mobilities are subjected to the same electric potential, they tend to move at 
different speeds and thus separate out . This is the basic principle of electrophoresis[18]. 
Electrophoresis, particularly in gels, affords a higher resolving power for biopolymers 
than any other method, and is therefore universally used for the analysis of protein 
preparations[19]. Figure 1.1.1 shows the electrophoretic pattern of blood serum 
proteins of normal and diseased people which is useful as a clinical tool. 

However, current methods in electrophoresis separation suffer from disadvantages 
like complicated design, construction, and operation of apparatus, low production 
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capacity, precipitation of proteins at high concentrations, etc. Evidently more research 
has to be done to make electrophoresis commercially viable. Hence simulation of this 
phenomenon can aid in further development . 

Specifically, the following goals can be accomplished : 

• Simulation can be used to correlate mobility/dissociation data with experimental data. 
A typical problem in this field is the general lack or inaccuracy of data. Using 
computational models can help to narrow-down inaccuracies . 

• Computational models reveal details of boundary structures like concentration, pH, 
and conductivity not easily amenable to direct experimental observation. 

® Simulations can be used to develop new techniques which are not yet available 
experimentally. This may be used to develop commercial high volume production 
techniques. 

® Numerical models can provide a rational basis for understanding many separation 
techniques. Existing theories are too restrictive in the sense that one theory can only be 
used to explain a particular technique. However, a mathematical model based on the 
governing equations can be simulated and can predict all techniques ( as will be 
illustrated in Chapter 4 ). 

Briefly, the thesis outline is as follows: 

Chapter 2 gives details on the development and testing of Petrov-Galerkin schemes. 
Chapter 3 outlines the algorithm used for flow field simulation along with test cases. 
Chapter 4 gives a detailed description of various electrophoretic techniques and their 
simulation. 

Chapter 5 gives the conclusions for Chapters 2, 3, and 4. 
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CHAPTER 2 

DEVELOPMENT OF PETROV-GALERKTN SCHEMES FOR 
CONVECTION-DIFFUSION PROBLEMS 


In this chapter a new Petrov-Galerkin finite-element formulation has been developed for 
transient convection-diffusion problems. Most Petrov-Galerkin formulations take into 
account the spatial discretization and the weighting functions so developed give 
satisfactory solutions for steady state problems. Though these schemes can be used for 
transient problems, there is scope for improvement. Basically the weighting functions 
used for the steady state problems are obtained by perturbing the functions which would 
otherwise lead to a Galerkin formulation. These perturbations are dependent only on the 
spatial discretization. However, for transient problems it was felt that this perturbation 
should also depend on the temporal discretization. The schemes proposed here, which 
take into account temporal as well as spatial discretization, provide improved solutions. 
In view of the generality of the differential equation being solved, these schemes can be 
used to solve any physical problem which is governed by the transient convection- 
diffusion equation. 


2.1 PROBLEM STATEMENT 


The transient convection-diffusion equation in n sd space dimensions is given as 

<(>, t + u • V<f> = V • ( K • V<> ) + / on Q (2.1.1) 

where § = <j) ( x,t ) is the dependent variable, a function of the spatial position x e Q 
c 9! nsd and time t e ] 0,T [ where T > 0 . The velocity field u = u ( x ) is given , 
and K is the conductivity matrix. The source term is given as / = / ( x ). The 
boundary T of the domain Q is assumed to have the following decomposition: 


r = r u r. 

8 h 


( 2 . 1 . 2 ) 
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0 = r s n r h ,2.1.3) 

where T and represent the Dirichlet and Neumann type boundaries, respectively. 
Hence, 


4> ( x,t ) = g ( x,t ) 


V x e T , t g ] 0,T [ 


(2.1.4) 


n( x ) • k • V <]>( x,t ) = /»(x,t) V X € r. , t G ] 0,T [ 


(2.1.5) 


where n is the unit normal vector to the boundary, and g and h are prescribed 
functions. 

The initial condition is 

- <{> ( x,0 ) - (J> 0 (x) x g Q (2.1.6) 

where <j> 0 is a given function.The initial/boundary-value problem is to find <{>=(()( x,t) 
which satisfies eqs. (2.1.1) and (2.1.4)-(2.1.6). 

2.2 FINITE-ELEMENT FORMULATION 

Consider a discretization of Q into element subdomains O e , e = 1,2,..., n el where n el is 
is the number of elements. Let r* denote the boundary of O e . We assume 

— n el 

Q = u Q e (2.2.1) 

e=l 

n el 

0 = n O e (2.2.2) 


The interior boundary is defined as 
r «- 2 n - r 


(2.2.3) 


Let V and S denote the finite-dimensional variational and trial-solution spaces such that 
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V = { wlwe H ! (Q) , w ( x ) = 0 V x e r } (2.2.4) 

& 

S = { O | O e H'(Q) , 0( x ) = g(\) V x e r } (2.2.5) 


It is also assumed that both the spaces consist of the typical C° finite-element 
interpolation functions. 

The discrete variational form of eq. ( 2.1.1) and the associated initial condition in eq. 
(2.1.6) are 

J {wO, t + wU*VO + Vw • (K • VO) } dQ 

+ f 1 J p{0„ + U*VO - V • (k • VO) -/ } dn 
e=1 n e 


= J w/zdr+ J w/d£2 (2.2.6) 

r h a 

J w(O-O 0 ) dn = 0 (2.2.7) 

a 

where p is a C* 1 perturbation to the weighting function w. The Euler-Lagrange equation 
corresponding to eq. (2.2.6) is 


if 1 J w{0, t + u*VO -- V©(k®VO) -/ } dQ 
e ~l Q e 

+ J w{ n*K*VO - h } dT + J w [ n®K®VO] dQ = 0 (2.2.8) 

^ h ^int 

where [ ] is the "jump" operator. 


Remarks 

(1) w is obtained by perturbing the weighting function w which when 


unperturbed leads to Galerkin formulation, i.e., the modified weighting function 
is given by 
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w = w + p (2.2.9) 

(2) The perturbation component weighs only on the element interiors and hence 
does not affect the boundary term h . 

(3) Brooks and Hughes [10] have shown that for rectangular elements the 
perturbation term does not affect the weighting of the diffusion term and it is 
expected that for reasonable element shapes the contribution is negligible. 

(4) The unknown (j) is interpolated with multilinear isoparametric interpolation 
functions. 

A consistent finite- element spatial discretization of equation (2.2.8) leads to the 
following set of ordinary differential equations: 

M d + C d = F (2.2.10) 

where M, C , and F are the mass matrix, stiffness matrix, and the generalized force 
vector respectively; and d and d represent the dependent variable and its temporal 
derivative at the nodes. 

This initial value problem is solved by a family of Predictor/Multi-conrector algorithms 
proposed by Hughes et al. [20]. 

The proposed weighting functions are the following: 

SUPG - a modified version of Streamline Upwind/Petrov-Galerkin formulation 
proposed by Brooks and Hughes[10]. 

TW - a scheme based on Transport Weighting proposed by Hughes et al.[12]. 

SW - a scheme based on Sigma Weighting proposed by Hughes et al.[12]. 

For all the above schemes the weighting function associated with element node 
number a is expressed as 

N a = N a + zP a (2.2.11) 

where N a is the multilinear weighting function (associated with node a) leading to 
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Galerkin formulation. 

The parameter z depends on the element Peclet number and is given [10] as 


[ a/3 0 < a < 3 

z = (2.2.5) 

l 1 a > 3 


where a e [ 0,°o ) is the element Peclet number( = uunh / 2k , where u is the flow 
velocity, h the element length along the advection direction, and k the diffusivity). This 
expression for z is a doubly asymptotic approximation of the form which leads to 
nodally exact solutions in one -dimensional steady state problems [10]. For convection- 


dominated flows z = 1 as a -»«*>. - 
The perturbation function P a for steady state cases is given as 


* 



where 


(2.2.13) 


P. = N* - N a (2.2.14) 

* 

Here N a is the weighting function for pure convection (i.e., for z = 1). The Streamline- 
Upwind/Petrov-Galerkin formulation given in [10] defines P a as 

P a = (h/2)s*VN a (2.2.15) 


where h is the element length along advection direction s. 

* 

The Transport Weighting formulation proposed in [12] defines N a as follows: 

N a ( x ) = Na ( x ) , x e d (2.2.16) 

x=x + es, e>0 (2.2.17) 

Here x represents the intersection with the element boundary of the line emanating from 
x in the direction s. 
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The Sigma Weighting formulation proposed in [12] defines N as follows: 


N a = cr /a , xe a 


a = < s • VN > 

a a 


<J = CT (2.2.: 

a = 1 a v 

where n en is the number of element nodes, Q e is the element domain, and < • > is the 
Macauley bracket ( i.e., <y> = y if y > 0, else <y> = 0 ). 

Note that P a for SW and TW is derived from eq.(2.2.14). 

For SUPG, P is given as 


P = x u«VN 

a a 


where x is a parameter having the dimension of time [10]. Based on x t an 
"algorithmic Courant number" [11] can be defined as 


C 2x = ( 2x ) I| U |1 


For transient problems, we propose 


P a = C 2t P a 


where C 2x , in general, is selected to be different from unity. For steady state problems, 
SUPG formulation has C 2x = 1. Hughes and Tezduyar[l 1] have used 

C 2t = C At (2.2.2^ 


for various problems governed by hyperbolic equation systems where C At is the element 
Courant number based on the time step, i.e.. 


C At = (At)|iu 


(2.2.25) 
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where At is the time step of the transient algorithm. In the next section a Neumann 
stability /accuracy analysis is carried out. This analysis is useful in arriving at estimates 
for the algorithmic Courant number c 2t . 

Remark 

Combining eqs. (2.2.11), (2.2.14), and (2.2.23), we obtain the expression, 

N a = (l-zC 2 ,)N a + zC 2T N a , (2.2.26) 

which is convenient for implementation purpose especially for weighting functions SW 
and TW because, unlike SUPG, for these weighting functions there is no explicit 
definition for P a . 

a 

2.3 STABILITY AND ACCURACY ANALYSIS 

Tezduyar and Hughes [21] have performed a detailed phase accuracy and damping ratio 
analysis for one-dimensional transient pure advection. In this analysis, the exact 
solutions were assumed to be of the form: 

<j>( x,t ) = e -^ + im>t e- ikx (2.3.1) 

and approximate linear finite-element solutions are of the form: 

-(£ + ioo)nAt -ikx« 

<KxA,tn) = e e (2.3.2) 

where % and \ represent exact and approximate algorithmic damping ratios (ADR) 

respectively; CD and co represent the exact and approximate frequency ratios(FR); 

and k denotes the wave number (k = 2ir fX where X is the wavelength). A dimensionless 
wave number can be defined based on k and the element length h as follows: 


q = k h 


(2.3.3) 
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The ADR is related to the amplitude decay of a wave while the FR indicates the presence 
of trailing or leading waves which are seen as oscillations. A good scheme should 
display minimum damping with a FR close to unity. Using these attributes as the 
guidelines and using the analysis in [21] , various plots were obtained for ADR vs q and 
FR vs q for a set of element Courant numbers, different weighting functions, and 

various forms of C 2x . The optimal value of for a given C A , was retained, and from 

this optimal data, C 2x was constructed as a function of C At . 

Figure 2.3.1 shows various weighting functions for one-dimensional linear elements 
using various forms of C 2x ( for z = 1 and C At = 0.5 ). 

Figures 2.3.2 to 2.3.7 show the ADR and FR characteristics of various forms of C 2x 
for SUPG and TW(or SW) weighting functions. 


Remarks 

£1) The various forms of C 2x are expressed as polynomial functions of C At , i.e., 

c 2t = a + ( 1 - a ) cf At 

where a and n are constants to be determined. 

(2) The forms C 2x = 1 and C 2t = C At have been extensively studied in [21]. 

(3) TW and SW weighting functions are identical in one-dimensional cases; 
for C 2t =1 they correspond to the Pade finite-difference approximation [22], 

Implicit Algorithms 

(1) Although algorithmic damping is zero for Galerkin formulation ( i.e., C 2x = 0 ), 

a rapid departure of FR ( not shown ) from unity predicts oscillations in 
numerical solutions. 

(2) For the SUPG weighting function, the best FR response is for the form, 

C 2t = 2/V15 + ( 1 - 2Wl5 ) C Al , 0 < C At < 1 (2.3.4) 
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with a reasonable ADR. As C At approaches zero, this scheme reduces to the 
one which has fourth-order phase accuracy [23] for the semi-discrete equations. 

(3) For the SUPG formulations, the forms C 2x = 2/Vl5 and C 2t = C At lead to 

satisfactory FR response. However, the corresponding ADR is not as good as 
for the case mentioned in remark (2). 

(4) For the TW( or SW ) weighting function of the form 

C 2t = 1/4 + ( 1 - 1/4 ) c£ t , 0 < C At < 1 (2.3.5) 

yields superior FR response. This form satisfies the unit CFL condition [24]. 

Explicit l-t>ass Algorithms 

(1) Galerkin ( not shown ) is unconditionally unstable. 

(2) SUPG ( C 2x = C At ) is second-order accurate. 

(3) SUPG ( C 2x = 1 ) and TW ( q t = 1 ) are identical. 

4 

(4) TW ( C 2t =1/4 + ( 1 - 1/4 ) C At ) is conditionally stable; the stability limit is 
C At = 0.253. 

Explicit 2-pass Algorithms 

(1) Galerkin ( not shown ) is unconditionally unstable. 

(2) All algorthms are second-order accurate. 

(3) TW ( C 2x =1/4 + (1-1/4) C A ^ ) is conditionally stable; the stability limit 
is C At = 0.303. 


2.4 NUMERICAL EXAMPLES 

In an effort to verify the improvements made in the weighting functions, a number of 
standard test problems were solved for convection dominated flows in one and two 
dimensions for various combinations of the weighting functions, algorithmic Courant 

number ( not to be confused with the element Courant number which is based on At), 
and the element Courant number. For one-dimensional problems even though the results 
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shown in figures are for C At = 0.6, the related discussions and conclusions also account 
for different value of C At . The problems are classified as follows: 

Advection of a Cosine-wave in One Dimension 

This consists of a cosine-wave initially on the extreme left and Neumann type boundary 
condition on the extreme right. The diffusion is set to zero, and the mesh and the 
convection velocity remain fixed. Tests were conducted for various time steps (i.e., for 
different element Courant numbers ). The exact solution consists of pure advection of the 
initial condition to the right. Solutions are compared at the end of the same time interval. 

Figure 2.4. 1 depicts the results for various weighting functions, various forms of C 2x , 

and a fixed value of C At . Since weighting functions TW and SW are identical in one 

dimension, only TW has been shown. Moreover, when C 2x — 0 , all weighting 
functions lead to Galerkin formulation. 

For SUPG for all forms of C 2v the peak amplitudes are comparable, and only trailing 

waves are observed. When C At is decreased, all forms of C 2x yield comparable 

solutions, and exhibit both leading and trailing waves. When C At is increased, the peak 

amplitudes are comparable; however, for the C 2x = 2/V 15 + ( 1 - 2/Vl5 ) C At form 
the trailing wave amplitude is halved compared to the other forms. For SUPG 

formulation, the form C 2x = 2/Vl5 + ( 1 - 2/V 15 ) C At performs better for high or low 

Courant numbers, while for medium Courant numbers, all forms of C 2x perform 
similarly. 

For TW (or SW) the form C 2x =1 produces leading oscillations; for other forms of 
C 2t the solutions are reasonable and comparable. A similiar trend is observed when the 

Courant number is decreased. With increasing Courant number all forms of C 2x give 
better results, and when C AI = 1, the solutions are identical. For TW(or SW) formulation 

_ 4 

the forms C 2x = C At and C 2x =1/4 + (1-1/4) C At provide satisfactory results for 
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a wide range of Courant numbers. 

Advection of a Discontinuity in One Dimension 

This is similiar to the previous problem of one-dimensional cosine-wave advection 
except that the initial profile has a discontinuity spread over one element. The exact 
solution again consists of pure advection of the initial profile to the right. Figure 2.4.2 

shows these results for C At = 0.6. For the SUPG formulation, all forms of C 2x provide 

comparable solutions (the peak overshoots are 7-9.5% of initial peak value and the 
undershoots are 0.7-2% of the initial peak value).When the Courant number is decreased 

the overshoots are higher for the C 2 T=C At form, and oscillations occur. This is not 

surprising because as C At tends to zero the scheme approaches Galerkin formulation. 
As the Courant number is increased, the solutions converge and become identical for 

C At = 1 for all forms of C 2r For the SUPG formulation all forms of C 2x yield 
reasonable solutions for medium to high Courant numbers. For low Courant number, 

the forms C 2x = 1 or C^ = 2/Vl5 + ( 1 - 2/Vl5 ) C Al perform better. 

For the TW( or SW) formulation, C 2x = 1 produces severe leading oscillations, while 

for the C^ = C At the undershoot is too high. However, much less overshoot and 

undershoot occur for the form C 2x =1/4 4 - (1-1/4) C At . For TW( or SW) the form 

4 

C 2x =1/4 + ( 1 - 1/4 ) C At performs optimally for the entire range of Courant 
numbers. 

Translation of a Cosine-hill in Two Dimensionsf Translating Puff ) 

The domain is a square region of size lxl in the positive Xj-x 2 plane with 30x30 

elements. A cosine-hill initial profile, distributed over 12 elements, is centered at 
coordinates ( 0.27, 0.50 ). The diffusion is set to 10' 6 . The velocity is unity in the 

Xj-direction and the time step is adjusted to give c* of approximately 0.25 at the tip of 

the hill. All boundaries are of Dirichlet type and homogeneous except at Xj =1.0 which is 

of Neumann type and homogeneous. The problem was run for 120 time steps, in which 
time the hill should be transported out of the region with minimal amplitude loss and 
minimal oscillations in front or in the wake. 
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Figures 2.4.3, 2.4.4, and 2.4.5 show the results for different weighting functions and 
various forms of C 2x . 

For SUPG formulation at step 80 it is observed that the peak amplitude is 99.2% ( as 
a percentage of the initial peak amplitude ) for C 2x - C At , 96.8% for C 2x - 1, and 

97.9% for the C 2x = 2/V 15 + ( 1 - 2/Vl5 ) C At form. The trailing wave amplitudes 
are comparable ( < 2% ). 

For the TW (or SW) formulations and all forms of C 2x the solutions are vastly 
superior to those obtained by using SUPG formulations( Peak height drop is less than 
1% while trailing wave amplitude is less than 0.01%). Notice, however, that for C 2x = 1 
the solution is marred by the presence of leading oscillations. 

For cosine-hill advection in two dimensions, the TW( or SW ) formulations perform 

_4 

quite well when the form C 2x = C At or C 2x =1/4 + (1-1/4) 0 At is used. 

Rigid Rotation of a Cosine-hill in Two Dimensions( Rotating Puff ) 

The setup conditions are the same as in the previous problem of two-dimensional 
cosine-hill advection except that the velocity field is rotational with respect to the center 

of the domain ( i.e., u : = -x 2 and u 2 = x l ), and all boundaries are of Dirichlet type and 
homogeneous. The problem was run for 200 time steps during which period the hill 
should perform one complete rigid body rotation. 

Figures 2.4.3, 2.4.6, 2.4.7, and 2.4.8 show the results for different weighting 

functions and various forms of C^. 

For the SUPG formulation the form C 2x = C At is apparently superior in terms of peak 

amplitude which drops less than 2% after 200 steps. For C 2x = 2/Vl5 + ( 1 - 2/Vl5 ) C A. 
form, the peak amplitude drop is less than 5% for the same trailing wave amplitude 
(<2%). 

For the SW formulation, the forms C 2x = C Ax and C 2x =1/4 + (1-1/4) C? At 
give comparable results( the peak amplitude drop is only 2% and trailing wave amplitude 
is less than 2% after 200 steps). However, C 2x = 1 fails totally for this problem, which 
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is not surprising in view of the one-dimensional problem results. 

For the TW formulation, the results are comparable to those of the SW formulation 

except that the form C 2t = 1 does not fail totally; however, the latter form produces 
leading waves with an amplitude of 9.5% which is quite unacceptable. 
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CHAPTER 3 

APPLICATION TO FLOW FIELD SIMULATION BY VORTICITY-STREAM 

FUNCTION APPROACH 

Flow field simulations for two-dimensional incompressible flows are generally based on 
either the primitive variable approach or on the vorticity- stream function approach. By 
using the vorticity- stream function formulation, the number of unknowns and transient 
equations are reduced, simplifying the problem for computational purposes. Since the 
governing equations involve transient convection- diffusion terms, the Petrov-Galerkin 
schemes developed in Chapter 2 were used for this formulation. 

3.1 VORTICITY-STREAM FUNCTION FORMULATION 

The fundamental equations for two-dimensional, incompressible flow of a Newtonian 
fluid with no body forces, and constant properties are the two momentum equations 
(Navier-Stokes) and the continuity equation [15], 

3u/3t + u 3u/3x + v 3u/3y = - 1/p 3P/3x + v ( 3 2 u/3x 2 + 3 2 u/3y 2 ) (3.1.1) 

3v/3t 4- u 3v/3x + v 3v/3y = - 1/p 3P/3y + v ( 3 2 v/3x 2 + 3 2 v/3y 2 ) (3.1.2) 

3u/3x + 3v/3y = 0 (3.1.3) 

These equations are in the primitive variable form involving velocity components, u, 

v and pressure P, and the fluid properties of mass density, p, and kinematic viscosity, v. 
The equations have been written in an Eulerian frame of reference and although it is 
possible to obtain numerical solutions from these equations[16], they can be further 
simplified both in terms of the number of dependent variables , and the number of 
transport equations being solved . This simplification leads to the vorticity-stream 
function formulation which has been used successfully for numerical simulations. By 
appropriate vector and algebraic operations on eqs.(3.1.1), (3.1.2), and (3.1.3), and 
introduction of vorticity and stream function as the new variables, the above three 
equations reduce to the following equations: 
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Vorticitv Transport Equation 

co, l + u*Vco = v V 2 co on Q. (3.1.4) 

where co is the vorticity, and u is the velocity vector. 

The boundary F is assumed to have the following decomposition , 


r=r ? u r r u r G 

0= r. n t s , 0= r- n r G> 0 =r G n r f 

Hence, 

co(x) = |'(x) V x e r~ 

v n( x )*Vcd( x ) = £(x ) V x e 


(3.1.5) 

(3.1.6) 

(3.1.7) 

(3.1.8) 


Nothing can be specified about the vorticity on the no-slip boundary, . Here n is the 
unit normal vector to the boundary, and g and h are prescribed functions. 

The initial condition is 

cd(x,0) — cd q (x) x e Q (3.1.9) 

Remark 

Vorticity is defined asVxu, a vector quantity. In two dimensions it reduces to a 
single component and can be treated as a scalar. 

Poisson's Equation for the Stream Function 

V 2x F + co = 0 on a (3.1.10) 

where is the stream function. 

The boundary T is assumed to have the following decomposition, 

r = q u r A u r G 

0= r, n r h , 0= r g n r G> 0 = r G n r h 


(3.1.11) 

(3.1.12) 



Hence, 
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¥(x) =s(x) 


V x e r 

g 

(3.1.13) 

n*V *F( x ) = h (x) 


v x e r. 

(3.1.14) 

*F( x ) = G (x ) 
n*V v F( x ) =H(x) 


V X 6 r c 

(3.1.15) 

) 




Remark 

The stream function *F is defined such that u = dW/dy and v = -B'F/Sx. Thus in 
eq.(3. 1.1 5), H is the wall velocity parallel to the boundary T c . 

It is required to find co = co(x,t) and V F = *F(x,t) for the two coupled partial differential 
equations. 

The flow is thus governed by the nonlinear vorticity equation, which is parabolic in 
time, and the elliptic Poisson equation for the stream function. For steady state solutions, 
which are of interest in many cases, the temporal derivative can be eliminated. However, 
interestingly, most numerical studies concerning steady state flows are based on the 
time-dependent equations, the steady state solution being obtained as an asymptotic limit 
of the unsteady state equation. This procedure allows the treatment of steady flows as 
parabolic in time, so that the solution marches in time to the steady state result. 
Frequently, the elliptic stream function equation is also treated as a time-dependent 
problem, with the steady state solution obtained for large time steps [25]. 

It is also worth noting that the vorticity transport equation serves as a model for many 
other transfer processes and the techniques used to solve it can be extended to application 
for processes including compressible flow. 

3.2 BOUNDARY CONDITIONS 

On no-slip boundaries both components of the velocity vector are specified. Hence it is 
necessary to specify boundary conditions for vorticity which are consistent with the 
stream function boundary conditions. Figure 3.2.1 shows the wall details at a 

boundary. If *F is expanded by a Taylor series out from a wall which is parallel to 

x-direction and where the values are denoted by 'F 1 , we get 
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= + )i Ay + l/2( 3 2x F/ 3y 2 ) x Ay 2 + l/60 3v F/ 3y 3 )j Ay 3 

+ 0[ (Ay) 4 ] (3.2.1) 


From the no-slip condition specified at the wall in eq.(3.1.15), ( dW/dy ) 2 = H , and 
9v/3x = 0 . 

Also 3u/3y = d 2x ¥/ dy 2 . Thus for a wall moving along x direction with velocity 
//, the vorticity is 

-co 1 = 3u/9y - 3v/3x = (2/h 2 ) ( *F 2 - X ¥ l -H h ) + 0[ h ] , (3.2.2) 

where h is the normal distance from the wall node to the adjacent node. Higher order 
forms have also been implemented though they tend to cause instability at high 
Reynolds number[26]. For stationary walls, eq.(3.2.2) reduces to 

co 1 = du/dy - 3v/9x = (2/h 2 ) ( X ¥ 1 - ^ ) + 0[ h ] (3.2.3) 

Remarks 

£1) On no-slip boundaries the boundary value of the vorticity depends on the stream 
function values and thus are time-dependent. 

(2) The above derivation is based on element boundaries normal to the no-slip 
boundary. This restricts the type of elements which can be used adjacent to the 
boundaries. The boundary conditions on vorticity can be imposed in a more 
consistent way within the finite-element context by a proper weak formulation of 
the problem. We were reminded of this by Glowinski[27]. 

(3) On outflow boundaries, a homogeneous Neumann boundary condition is 

usually specified for co and T / . Such a "traction-free" boundary condition, 
although suitable for Navier-Stokes equations, is not the best one for the 
vorticity-stream function formulation. A more appropriate "absorbing" boundary 
condition has been proposed by Glowinski[28]. 
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33 FINITE-ELEMENT FORMULATION 


Consider a discretization of £2 into element subdomains £2 e , e = 1,2,..., n el , 
where n el is the number of elements. Let T 6 denote the boundary of Q e . We assume 

n = v a e (3.3.i) 


n el 

0 = n n e 

e=l 


(33.2) 


The interior boundary is defined as 

r hu = w n - r (3.3.3) 

e=i 

Weak form of the Vorticitv Transport Equation 

Let V and S denote the finite-dimensional variational and trial solution spaces such that 

V = { w I w € H ! (£2) , w( x ) = 0 V x e , T G } (33.4) 

S = {co|co e H 1 (Q),co(x)= g(x) V x e T-, 

co(x)= C0 G (x)V xe r G ) (33.5) 

where C0 G represents the value attained by vorticity on boundary V G . 

The variational form of eq. (3.1.4) and the associated initial condition in eq.(3.1.9) 
are 

| {wco, t + wu*Vw+ v Vw*Vco) d£2 + 

Q 

n ei r 2 r 

£ J p {co, + u*Vo) - v Vco}d£2 = J w/TdT dT 

Vw e V (3.3.6) 

J w (co-fflo) da = 0 V w e V (3.3.7) 

Q 
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where p is a C' 1 perturbation to the weighing function w, i.e., the modified weighting 
function is given by 

w = w + p (3.3.8) 

The Euler-Lagrange form of the above equation can be derived and is similiar to the 
one derived in Chapter 2, i.e.. 

If 1 J w{co„ + u»Vco - vV 2 0) } dn 

e=l n e 

+ J w{ v n»Vco- h } dF + J w [ v n®Vco] dQ =0 (3.3.9) 

T h r int 

where [ ] is the ’'jump" operator. 

Weak Form of the Poisson’s Equation 

Let V and S denote the finite dimensional variational and trial solution spaces such that 

v = (w i w g h ! (D) , w( x ) = o v x g r , r G } 

S = {^|^ G h ! (Q) , 'f ( x ) = # ( x ) v x G r , 

y ( X ) = G ( x ) V XG r G } 

The variational form of eq. (3.1.10) is 

J Vw»V x F df2 = j woo dH + j w/z dT V w g V (3.3.12) 

a q. r h 

where w is the weighting function leading to Galerkin formulation. 


(3.3.10) 

(3.3.11) 
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3.4 MATRIX FORM OF THE DISCRETIZED EQUATIONS AND THE 

ALGORITHM 

A consistent finite-element spatial discretization of eqs. (3.3.8) and (3.3.12) leads 
to the following set of ordinary differential equations: 

M d + S d = F (3.4.1) 

Kf = F • Md (3.4.2) 

where M and M are the mass matrices, fc and K stiffness matrices, F and F the 
generalized force vectors, d and d the vorticity and its temporal derivative at the nodes, 
and the stream function at the nodes. Note that the mass matrices in eqs. (3.4. 1) and 
(3.4.2) are not identical since the weighting functions used for the equations come from 
different spaces. 

Predictor/Multi-Corrector Algorithm 

The algorithm used to solve eqs. (3.4.1) and (3.4.2) is a modified version of the 
Predictor/Multi-corrector algorithm used by Hughes and Brooks[20] and is shown 
below: 


(1) Predictors 
( 0 ) 

d n + l = d „ + At(l-Y) d n 

• ( 0 ) 

<W = 0 

Set i = 0 and go to step (2a) or (2b) 

Predict the vorticity at next time step based on the current step. 


(2a) Solution alternative I 


(0 e e (i) e 

CY n + i = A { f - m (d n+1 ) } 

0) (i) 

( d c >n*l = L 


Stream function 
calculation 

Update vorticity boundary 
condition 
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(5) Convergence check 

If the vorticity equation residue norm is within the tolerable limit, then increment 
n by 1 and go to step (1). Else increment i by 1 and go to step (2a) or (2b). 

If the vorticity has converged, then go to next time step else keep correcting it. 


Remarks 


(1) At is the time step chosen and y is the Euler parameter chosen between 0 and 1, 
while subscript n represents time levels and superscript i represents the iteration 

level. Symbol A represents assembly over the elements e= l,2,..n el . 

(2) M* is defined as follows: 

M = M + At y K for implicit solutions 
M* = M l for explicit solutions 

where M L is the lumped mass matrix. 

(3) For step (2a), if N is not large enough the solution diverges. From practical 
experience it is seen that convergence is slow with step (2a) and hence step (2b) 
is used . 

(4) Weighting functions developed in Chapter 2 are used for the vorticity transport 
equation since they have good stability/accuracy characteristics. 

(5) Matrix L is obtained from eq. (3.2.3) which relates the wall yorticity( d G ) with 
the stream function. 

(6) Since matrices M and K for the vorticity equations depend on the spatial and 
temporal discretization besides the flow field they are reformed every iteration 
within each time step. However matrices M and K for the stream function are 
formed only once in the first iteration of the first time step since the weighting 
function depends only on the spatial discretization. 
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3.5 NUMERICAL TEST CASES 


Flow Past Step 

This is a relatively stringent problem for testing the algorithm because of sharp 
boundaries which result in steep gradients of variables. Convergence in this problem is 
thus not expected to be as good as for problems involving smooth geometries like flow 
past a cylinder. ( The latter case and associated phenomenon of vortex shedding is 
currently being studied by Glaisner[29], and, results obtained are in strong agreement 
with those obtained by Fomberg[30], Smith and Brebbia[31], and Brooks and 
Hughes[32]). Figure 3.5.1 gives a problem description and the appropriate boundary 
conditions used. This problem was solved using a crude(210 elements) mesh and a fine 
mesh (804 elements) for Re = 200 . The mesh was not refined near the walls in order to 
study the effects of a uniform mesh. 

The solutions obtained for Re = 200 cases are shown in Figures 3.5.2 and 3.5.3 . 
Results obtained compare quite well with those obtained by Hughes et al.[32] using a 
penalty formulation. However, in our computations the Dirichlet boundary values are 
allowed to reach the full values after 10 time steps. Thus solutions obtained should be 
compared with caution. 

The following observations can be made: 

(a) As the mesh is refined, the time step has to be reduced even for the implicit 
calculation because of the evolving vorticity boundary condition on the walls, which is 
proportional to 1/h 2 . This is a limitation because an accurate solution obtained by a fine 
mesh would take more time steps to compute. 

(b) In general, the viscous effects dominate at the walls resulting in high wall vorticity. 
Away from the walls the flow dominates and convects away the vorticity. This suggests 
mesh refinement near the walls. 

(c) For Re = 200, a fine mesh with low time-step results in a region, of recirculation 
downstream of the step which develops further in time. This is in agreement with the 
solution obtained by Hughes, et al.[32]. 

Driven Cavity Problem 

This, like the step problem is a stringent test for the algorithm and is described in Figure 

3.5.4 with the appropriate boundary conditions. A uniform mesh of 20x20 elements was 
chosen and the solutions were obtained for Re = 100, 400, and 800. In all the cases, the 
time step was fixed at 0.2 and the problem was run for 40 time steps. Transient results 
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are shown for Re = 100 in Figure 3.5.5, andthe last step results for Re = 400 , and 800 , 
in Figures 3.5.6 and 3.5.7. The solution can he compared with Hughes et al.[32) who 

used 20x21 elements . Figure 3.5.8 shows the results obtained by Hughes et al.[32]. 
The following observations can be made: 

(a) Initially, the vorticity is generated by the top moving wall and is transported to the 
adjacent walls by the flow developed. However, as a consequence of this generation and 
transportation of vorticity, a circulating flow develops in the cavity. At steady state, the 
vorticity diffusion is balanced by the vorticity convected and recirculated. 

(b) As a result of the circulating flow in the cavity, two small regions of recirculation 
develop on the lower comers.These grow with Reynolds number until Re = 800. 

(c) In general, agreement is good for low Reynolds number with solution obtained by 
Hughes[32]. Note that both vorticity and stream function are practically constant after 
40 steps for low Reynolds number. 
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CHAPTER 4 

APPLICATION TO SIMULATION OF ELECTROPHORESIS SEPARATION 

PHENOMENA 


Before actual simulation of any process it becomes necessary to identify the governing 
equations (PDE’s, ODE’s, algebraic equations, etc.) based on certain basic laws (mass, 
momentum, energy, charge balance, etc.) and also to understand the 
limitations/assumptions made in arriving at the equations. Electrophoresis simulation 
also requires the following: 

® Development of special techniques to solve problems in convection dominated 
flows. Successful separation of particles is based on large differences in mobilities and 
low diffusion rates of individual species. Evidently, as explained in Chapter 1, 
convection dominated flows need Petrov-Galerkin formulations which have been 
developed in Chapter 2. Thus the schemes developed in Chapter 2 have been used in 
electrophoresis simulation. 

• Prior knowledge of the external flow field (i.e., the velocity of the solvent in 
which the charged particles have been dissolved). In Chapter 3, a method has been 
developed to model the flow field based on the vorticity -stream function 
formulation. This formulation can be used to solve for the external flow field. 
However, it must be noted that in most electrophoretic separations the external velocity is 
kept at zero to avoid remixing of the separating species. 

4.1 CLASSIFICATION OF SEPARATION TECHNIQUES 

Although there are several classification of separation techniques the three distinct 
categories studied by us are listed below: 

zone electrophoresis ( ZE ) 

moving boundary electrophoresis ( MBE ) 

isotachophoresis ( ITP ) 


The above methods differ in the applied initial condition which leads to a different 
separation mode. 
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In £JL the whole separation column is filled with a single basic electrolyte system with 
a specific conductivity. The mixture of substances to be separated is injected at a certain 
location in the column and the passage of current mobilizes the substances at different 
rates and thus zones are created. Figure 4.1.1 illustrates this process schematically. The 
shape of the individual zones are influenced by the potential gradient If the electrolyte is 
such that its conductivity is not affected by the substances to be separated, then the 
potential gradient remains constant and it is the most common type of ZE ( this is the 
case simulated). In some cases the conductivity does change and the potential gradient is 
not a constant. Thus, when the species enters a region of high potential gradient 
(absolute ) it speeds up, leading to a spreading of the front boundary . When it exits , the 
boundary sharpens. The opposite is true for the back boundary. A similiar deduction can 
be drawn for species entering or exiting low potential gradient( absolute ) regions. 

Some characteristic features of ZE are the following: .. - 

(i) particles with either sign can be separated. 

(ii) concentrations decrease permanently. 

(iii) zone boundaries can sharpen, or spread. 

(iv) zone velocities are different for different species. 


MBE resembles ZE closely and is characterized by separation of the column into three 
parts. Both sides are filled with the same basic electrolyte; the middle contains the 
mixture to be separated. Since the middle width is comparable to the column length, 
complete separation into individual zones cannot be accomplished. In yet another 
variation of this method, the compartment may be divided into only two regions - one 
containing the basic electrolyte and the other containing both the electrolyte and the 
mixture to be separated. Figure 4.1.2 shows the latter mode schematically. In the 
simulation, both species A and B are positively charged with A having a higher mobility. 
Therefore the leading zone contains only A followed by a zone containing both A and B. 
Only species A can be fully separated. 

Some characteristic features of this method are the following: 

(i) practically used to separate particles of one sign. 

(ii) only one species is fully separated, others are mixed. 

(iii) boundaries may sharpen or spread. 

(iv) zone widths increase or decrease permanently. 

(v) zone velocities are different for different substances but remain constant. 
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TTP is characterized by the fact that a cation/anion mixture cannot be separated by this 
method. The separation column is divided into three unequal parts. One part is filled with 
a leader electrolyte, followed by the mixture to be separated, and finally, the terminating 
electrolyte. If the mixture has anions to be separated, then the leader anions must have 
the highest mobility and the terminator the lowest in the column. On passing current, the 
mixture separates completely into individual substances with the anions arranged in the 
order of mobilities. The leader is in the front followed by the species to be separated, and 
finally the terminator. After steady state has been achieved, all the boundaries move with 
the same constant speed . Figure 4.1.3 depicts this method schematically. 

Some characteristic traits of this method are the following: 

(i) complete separation of species is achieved after steady state is reached. 

(ii) concentrations can increase or decrease during transient phase. 

(iii) concentrations are constant after steady state achievement. 

(iv) zone widths increase or decrease during transient phase. 

(v) zone widths are constant after steady state achievement. 

(vi) zone velocities are same and constant after steady state achievement. 

The above 3 separation techniques are used in one-dimensional separation. In two 
dimensions, there is a different classification, 

zone electrophoresis - similiar to one-dimensional ZH 

cross electrophoresis - oppositely charged species are made to cross each other so 
that a chemical interaction produces a new species. 

diagonal electrophoresis - basically ZE along first direction followed by a chemical 
treatment of separated species, followed by ZE along second direction, 
cross flow-through electrophoresis - an external bulk velocity is applied along 
one direction and the electric current perpendicular to it. Separation is achieved 
continuously along each electrode. 

Remark Although there are several other criteria used to classify separation processes 
(amount separated, high/low voltage, continuous/batch etc.), the above classification 
methods are more amenable to mathematical modelling since they provide definite 
information regarding domain, initial concentration profiles, etc. 
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4.2 ASSUMPTIONS 

The governing equations are derived with the following assumptions: 

Fully ionized species - while charged or uncharged species can be easily handled by the 
model, the charged species are usually in equilibrium with their neutral counterpart. This 
represents a chemical reaction which complicates the model and is currently not being 
studied. 

Individual species do not affect each other - the basic assumption made in deriving the 
flux equations is that flux of one species does not depend on the other. This is accurate 
only in dilute solutions. 

Incompressible flow - since, in electrophoresis, one is dealing with common liquids at 
low velocities, this assumption is valid. However, if a chemical reaction changes the 
densities spatially and temporally, this assumption would have to be scrutinized. 

Electrical heating effects negligible - the passage of an electric current induces heating of 
the liquid column which can in turn, induce thermal convection currents and remix the 
separating species. For this reason the separation column cross-section is kept low. 

Electro-osmotic flow negligible - this flow is generated when an electric field is applied 
to the system and the charge in the diffuse part of double layer of the wall moves to the 
oppositely charged electrode. This drags particles close to the double layer along with it 
thus producing a net flow. 

Electroneutrality - i.e., the sum total of all charge at a point (at the continuum level ) is 
zero. Most fluids satisfy this property except in a region very close (100 Angstroms) to 
the electrodes. 



32 


4.3 GOVERNING EQUATIONS 

The flux of each dissolved species is given by[33] 

N^-z^Fc-VO - DjVcj + (4.3.1) 

electromigratory diffusive convective 

flux flux flux 

where 

Cj = unknown concentration of the i* species 

O = unknown electric potential 
z i = signed valence of the i* species 
u i = mobility of the i* species 

F = Faraday's constant 
v = bulk fluid velocity 
D i = diffusion coefficient for the I th species 
Rj = rate of generation of the i* species 
A material balance yields 

3c, 

— = -V®N, + R. ( 4.3.2 ) 

3t 

Substituting eqs. ( 4.3.1 ) in ( 4.3.2 ) and using the incompressible flow assumption 
( V®v = 0 ), we get 

+ v®V c { = ZjUjFV^CjVO ) + V^DjVc.) + Rj (4.3.3) 

3t 

or on simplification, 

2 .2 

+ ( v - z^F VO^Vcj = z^F V <E> + Dj V c { + R i 


3 Ci 

at 


(4.3.4) 
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Moreover the current density is given by 


= FI 


Z iN; 


(4.3.5) 


On multiplying eq.(4.3. 2 ) with ZjF and summing over all species, we get 


d 


-FSz i c i = -V*FSz i N 1 + FSzjR. 
9t 1 i i 


(4.3.6) 


Using the principle of electroneutrality and charge conservation ( I z i R i =0 and 

i 

I Zj Cj = 0 respectively ) gives 


i = -F 2 VO I z i 2 u i c i - FI z i D i Vc i 


(4.3.7) 


V • i = 0 

Further assuming a constant current density, we see that 
-VO = i/K + F/k I z i D i Vcj 


(4.3.8) 


(4.3.9) 


If electroneutrality is imposed explicitly, we get 
- V <D = i/K + F/k I Zj(D f D n ) V Ci 

ten 


(4.3.10) 


where 

K = F 2 Iz i 2 u i c i 
i 

If electroneutrality is imposed explicitly, we have 


(4.3.11) 


K = F 2 I ZjC; (ZiU f Z n U n ) 


(4.3.12) 
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Thus the equations to be solved are the following : 

+ ( v - z-u-F V<3>)*Vc- = z.u-Fc.V^ + D- V 2 c- + R- (4.3.13) 

3t 

where the electric potential is evaluated from 


-Vd> = i/K + F/K X Zi (D f D n ) VCj (4.3.14) 


and the conductivity from 

k = F 2 X z i c i (z i u r z n u n ) (4.3.15) 

i#n 


Remarks 

The above eqs. (4.3.13), (4.3.41), and (4.3.15) are thus used to evaluate the 
concentrations of n-1 species and the n* species concentration can be evaluated from 
electroneutrality, i.e., 

Vn = £ (4.3.16) 

i*n 


4.4 STRONG FORM OF THE PROBLEM 

The strong form of the problem based on the equations derived in previous section can 
now be written as follows ( similar to the ones derived in the previous chapters except 
that there are n unknown species concentrations to be evaluated at a given node): 

3c i 

+ ( V - Z-U-F VO)«VC: = ZjU:FC:V 2 <E) + D- V 2 C: + R- 

* 1 1 111 11 1 

3t 


on Q and i = l,2,..,n-l 


(4.4.1) 
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where 


- VO = i/K + F/k I z i (D f D n ) VCj 

i^n 

and 

K = F 2 X ZftCZjUj-Z^) 


(4.4.2) 


(4.4.3) 


The boundary T of the domain Q , is assumed to have the following decomposition, 


r = r g i u r A i (4.4.4) 

0 = n r^i (4.4.5) 


where F^i and T^i represent the Dirichlet and Neumann type boundaries, respectively. 
Hence, 


Cj( x,t ) = g x,t ) 


Vxe r i.t g ] 0,T [ and i = 1,2, ...n-1 (4.4.6) 
& 


n(x)*Vc i (x,t) = /i i (x ) t) Vxe T^i.t g ] 0,T [ and i = 1,2, ...n-1 (4.4.7) 


where n is the unit normal vector to the boundary, and g . and h . are prescribed 
functions. 

The initial condition is 

c i ( x,0 ) = c i0 ( x ) x g £2 , i = 1,2,. ..n-1 (4.4.8) 

where c i0 is a given function.The initial/boundary-value problem is to find Cj = Cj( x,t) 

which satisfies eqs. (4.4.1), and (4.4.6)-(4.4.8). It is assumed that the current density 
is known, i.e., i = i ( x,t ) is specified. 
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4.5 FINITE-ELEMENT FORMULATION AND THE ALGORITHM 

Adopting a procedure similiar to one in Section 2.2, a consistent finite-element 
spatial discretization of eq(4.4.1) leads to the following set of ordinary differential 
equations: 

Md + Kd = F (4.5.1) 

where all the symbols have their usual meaning. 

However solving this set of ODE's is different because the convection velocity induced 
due to the electric field has to be evaluated. In other words, these ODE’s can be solved 
only if the potential gradient is known at each time step. 

Thus the scheme used to solve this system is similiar to the Predictor/Multi-corrector 
algorithm used by Hughes and Brooks[20], except that several corrections are required 
to evaluate concentrations and the potential gradient. 

The potential gradient is evaluated at nodes based on a 2 nd order finite-difference 
scheme, i.e., in evaluating VO from the expression, 

-VO = i/K + F/k S Zj(D r D n ) VCj (4.5.2) 

a 2 nd order finite-difference scheme based on nodal concentrations is used to 
approximate the VCj term. 

Remarks 

(1) Since a finite-difference scheme is used to evaluate derivatives based on nodal 
concentrations, it restricts the usage to elements with boundaries parallel to the 
X and Y-axis. 

(2) Depending on the spatial position and time, the flow for a particular species 
might be convection dominated and for another species diffusion dominated. 

This situation is possible because the convection term for each species is 
governed by the potential gradient (common to all species and a function of 
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space and time), and its particular valence and mobility. 

(3) A Courant number based on "electro-migration*' can be defined which will be useful 
in stability/accuracy analysis, 


(EC At )j = (Fz iUi V<D)At 
h 


i = (4.5.3) 


Here (EC^ is the "Electro-migrative Courant Number" for the i* species based 
on time step At, and element length h, along electromigation direction. 


4.6 NUMERICAL TEST CASES 

Although it is not possible to predict beforehand the course of separation process, in 
1897 Kohlrausch established a general principle for understanding the general features 
which can be stated as follows: 

" The passage of an electric current through an electrolyte system causes 
changes only where the system is nonhomogeneous." 

Consequently, 

(i) the passage of electric current through a homogeneous solution leaves the 
concentrations unchanged. 

(ii) the passage of an electricxurrent causes the boundary between two ionic species to 
advance and either degrade or self-stabilize in time. 

(iii) concentration gradients are largely unaffected, but flatten out as a result of diffusion. 

A review of numerical simulations indicates a paucity of work done in the field of 
electrophoresis. The models developed in [34] and [35] are specialized for either MBE or 

ITP in one-dimensional steady state cases. The general approach to electrophoresis 
presentedby Bier et al.[36] does explain all methods in one dimension. In two 
dimensions, surprisingly, no work seems to have been done - this despite the fact that 
two-dimensional electrophoresis is a dominant method in resolving proteins. The 
algorithm developed by us has the following advantages: 

(i) has good stability and accuracy characteristics as demonstrated in Chapter 2. 

(ii) can handle several species. 
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(iii) can handle external flow fields. 

(iv) can handle irregular domains (however, see remark in Section 4.5). 

One-dimensional simulations Based on the classification given in Section 4.1, 
one-dimensional simulations have been done for ZE, MBE, and 3TP. The results have 
been compared with those of Bier et al.[36], and Fidler et al.[34]. However, Bier's 
simulations take into account chemical dissociation equilibrium which can affect the 
solution obtained depending on the dilution. 

In all the simulations the two species to be separated are denoted by A and B, and the 
basic electrolyte in which they are introduced consists of P + ( or L + and T 4 * ) and Q' ions. 

Zone Electrophoresis in One Dimension 

The initial condition consists of a uniform concentration distribution of P + (and Q‘) ions 
and a cosine-hill for both A + and B“ bearing opposite charges. 

On the application of an electric field, each species migrates a distance proportional to 
its mobility and charge. In this way the species separate. The problem was run for 150 

time steps ( At =0.02 ).The temporal evolution of A + ,B' ,P + , and potential gradient (V<|)) 
are shown in Figure 4.6.1. Broadening of the leading boundaries of A + and B~ indicates 

that the absolute value of V(j) is high off center .This is consistent with the V(|> plots. 

Since the potential gradient is a function of the current density, concentrations and 
concentration gradients, adding the mixture S will disturb the existing uniform potential 
of the electrolyte. This is seen as a cosine-hill hump in the potential gradient initially. 
However as the species separate, it attains a much more complicated profile than is 
shown for this simple case. 

On comparing this with the results obtained by Bier et al.[36] there is a similar trend 
observed except that the species to be separated has a sharpened leading boundary and 
smeared trailing boundary. This is because in [36] the species is an ampholyte ( i.e., its 
mobility is a function of other species) and the electromigration velocity is thus 

dependent on the V<J), as well as the variable mobility. This was also simulated by 
making the mobility of A + , a variable dependent on P + ion concentration in our 
computation ( not shown ) . 

Moving Boundary Electrophoresis in One Dimension 

The initial condition were chosen to compare the result with those of Bier et al.[36]. This 
initial condition is not the same as shown in Figure 4.1.2. The problem was run for 250 
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time steps and since B + has a higher mobility than A + , the latter separates out fully as 
expected. This is to be contrasted with the example shown in Figure 4.1.2 where species 
A + with highest mobility separates out fully. This indicates that the initial condition and 
particle charge play a significant role in the fate of separation. Moreover, the boundaries 
can shaipen ( self stabilize) or deteriorate ( diffuse ), as shown for the two boundaries. 

The temporal evolution of A + , B + , P+, and V<$) are shown in Figure 4.6.2 at regular 
intervals of 25 steps. Species A + is not affected too much by the current because it has 
a low mobility and hence, tends to diffuse. Species B + on the other hand, shows clearly 
the effects of electromigration with some initial profile adjustments. P+ shows a tendency 
to reach spatial peaks which spread out in time. The potential gradient undergoes 
changes which are difficult to predict but shows plateaus in regions where concentrations 
are more or less constant. 

These results are in general agreement with those of Bier et al.[36] for the descending 
arm case. However the comparison is qualitative, based on profiles. 

Isotachophoresis in One Dimension 

The initial condition for ITP is depicted in Figure 4.6.3. The mixture is introduced 
between the leader L + and terminator T* such that 
u L > u A .u B > u T 

Also note that the terminator T 1 " is not strictly separated from the mixture which is not 
in general compliance with the requirements of ITP mentioned in Section 4.1 . 

On passing current the species separate out in the order of their mobilities and finally, 
the separated species move as a single system . This problem was run for 1750 time 
steps and an intermediate and final condition are shown in Figure 4.6.3. Except for the 
leader, all other species undergo significant changes in zone width and peak 
concentrations. Moreover as a result of numerical smearing, the boundaries tend to 
overlap. 

The temporal evolution of A+, B+, T+, L+, and V<j) are shown in Figure 4.6.4 at 
regular intervals of 150 steps. ForT+, in the initial stages the concentration levels off to 
zero. However later it levels off to a positive value signifying approach of a steady state. 

Species B + follows a complicated evolution with increases, decreases, and plateaus in 
the concentration levels until it flattens out at steady state. 

Species A + follows a complicated evolution with increases, decreases, and plateaus in 
the concentration levels until it flattens out at steady state. However the steady state 
concentration level is higher than B + while the width is less. 
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The leader L + shows relatively little change except for an initial boundary profile 
adjustment. 

The potential gradient follows an initial period of adjustment until the species separate. 
After steady state is attained it develops local plateaus which correspond to the plateaus 
of each species. 

On comparing this simulation with that of Zidler et al.[34] the following differences 
emerge: 

(i) The initial concentration levels of A + and B + are not equal in [34]. 

(ii) Initially, the terminator T + is strictly separate from the mixture in [34]. 

The steady state results as obtained by us are in good agreement with those of Zidler et 
al.[34] except for the behavior of the terminator which is different because of the above 
two reasons. Moreover, the results by Bier et al.[ 36] also indicate that both increasing 
and decreasing profiles of the terminator are possible depending on the initial 
concentrations which is consistent with our results. 

Two-dimensional simulations As pointed out earlier in this section, there is no result 
available to compare with in two dimensions. However, in order to test the algorithm 
and to demonstrate its capability, results for zone and cross (no reaction) electrophoresis 
have been shown. Although other methods have also been simulated, the results have 

not been shown here. For both cases the domain was taken as 0.66 x 1.0 units in the 

positive XY plane with a mesh of 20 x 30 square elements. For simplicity, the constant 
current density was applied in the Y-direction only. 


Zone electrophoresis in Two Dimensions 

Species A + and B" , bearing opposite charges, are introduced as cosine hills centered at 
coordinates (0.33,0.5). The basic electrolyte, as usual, is made of P + and Q' ion 
distributed uniformly ( not shown ). 

The problem was run for 100 time steps and results plotted at steps 0, 50 and 100, in 
Figure 4.6.5. The profiles along Y-direction are similiar to those attained in the 
one-dimensional version of ZE. However, both diffusion and electro-migration act in 
the X-direction. Though diffusion is understandable, electro-migration might appear 
peculiar because, as stipulated earlier, the current was applied only in the Y-direction. 
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However, on careful examination of eq. (4.3.11) it becomes obvious that the potential 
gradient also depends on the concentration gradients which do exist along X-direction. 

This is confirmed by the nonzero d§/dx which develops in time (not shown). 

Cross Electrophoresis in Two Dimensions 

Species A + is introduced as a cosine-hill centered at coordinates ( 0.33,0.25) and species 
B' centered at ( 0.33,0.75) so that they cross over in the domain on applying current. 
This problem was run for 250 time steps and the results are shown at steps 0, 50, and 
200 in Figure 4.6.6. 

Species A + behaves quite normally, displaying the usual broadening and reduction in 
concentration level. However B + shows an increase in concentration which is accounted 
by the sharp potential gradient change at step 50( It does not affect A + because it has not 
reached there ). However this sharpness diffuses later which is consistent with 
Kauhlraucsh's principle. Moreover, it is apparent that the potential gradient is 
predominantly controlled by B + because it is seen to drift along with the profile of B + 
rather than that of A + . 
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CHAPTER 5 
CONCLUSIONS 


The results from Chapter 2 show that by choosing a proper weighting function and an 
appropriate form for the algorithmic Courant number, a variety of problems can be 

solved giving improved results. Choosing the form C^ — 2/Vl5 + ( 1 - 2/Vl5 ) C At 

with SUPG and C 2t = 1/4 + ( 1 - 1/4 ) C At form with TW (or SW) gives optimal 
results for a range of Courant numbers, for one-dimensional problems. For 

two-dimensional problems, in addition to the above forms of C 2x , the form C 2x = C At 

also performs well. In general, the extension of one-dimensional schemes to 
multi-dimensions should be done with care. However, the proposed weighting functions 
perform satisfactorily in multi-dimensions as is clear from the results. 

Based on the results obtained, and from the stability and accuracy analysis performed, 
one can deduce that the proposed transient schemes are versatile, produce wiggle-free 
solutions, and show minimal damping and a good frequency ratio response. 

Numerical test cases run for flow field simulations in Chapter 3, indicate that the 
vorticity-stream function algorithm behaves quite well at low Reynolds number. 
However, there is limitation on the time step, which depends on the mesh size and is 
thus crucial for the stability/accuracy of the solution. 

As Reynolds number increases, sharp changes occur in the solution which lead to 
resolution problems. Though these can be resolved by mesh refinement (which is an 
extra computational burden), in complex domains it is not possible initially to predict 
regions of steep gradients. 

For high Reynolds number this model is not appropriate because turbulent flows are 
three-dimensional in nature[37]. Also, resolving the small scale structures would require 
refinement beyond feasible limits [37]. 

Moreover, if the pressure is to be determined, the primitive approach, in terms of the 
velocity and the pressure, is often advantageous. For flow in finite-sized enclosures and 
flow over arbitrary three-dimensional bodies, this assumption of two-dimensional flow 
may not be applied and the general three-dimensional circumstance is to be considered. 
The problem is extremely complicated in that case, and though the vorticity-stream 
function approach may be extended to this flow, the primitive variable method is often 
more appropriate.Thus other approaches like penalty formulation, primitive variable 
approach, etc., should also be tested with the newly developed Petrov-Galerkin 
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schemes. 

The numerical scheme proposed for electrophoresis simulation is capable of simulating 
ZE, MBE, and ITP in one dimension. The results obtained are in good agreement with 
those obtained by others. In general, electric current application mobilizes the particles, 
which in turn, alter the potential gradient distribution. Thus the electromigration velocity 
will change as a result of the changing potential gradient. In many cases, complicated 
concentration and potential gradient profiles are achieved which are difficult to predict 
otherwise. 

However, a combination of Petrov-Galerkin finite-element and 2 nd order finite- 
difference discretization was used in the algorithm which does place a restriction on the 
domain shape. In order to avoid such a restriction, the algorithm should be modified for 
quadratic elements and thus a pure finite-element scheme can be used on any domain. 
This would require testing of the Petrov-Galerkin weighting functions which have been 
developed and tested only for linear elements. However, it is expected that the 
performance would be more or less the same. 

The effects of thermal heating/osmotic flow have been neglected in our model. By 
including such effects and learning to control them, it is possible to keep them to a 
minimum and thus develop better commercial techniques. 

Reaction terms have not been included as yet in the test cases; clearly, to take into 
account equilibria between neutral and charged species, this is essential. This would 
require a better understanding of individual physical problems. 

The assumption of a constant current density, while valid for one-dimensional 
problems, is not a good approximation for two-dimensional cases; it is expected that 
better models will be developed. Two-dimensional simulations are quite expensive but 
provide invaluable help in understanding the physical process. In short, more work 
needs to be done to simulate electrophoresis separation fully, but in the meantime the 
existing model does provide useful insight into the physics of the problem. 

To sum up, the motivation for. developing accurate time-dependent Petrov-Galerkin 
schemes was to simulate flow field and electrophoresis separation phenomena; 
nevertheless, the equations considered might be those governing thermal convection/ 
diffusion, chemical convection/diffusion/ reaction, etc. Our results indicate a strong 
potential for future work in this area and it is expected that such formulations will be 
used to obtain better transient solutions for Navier-Stokes and Euler equations which 
have become "classical'' challenges over the years. 
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a 


Figure 1.1.1 Electrophoresis. Characteristic examples of separation of pathological blood 
serum proteins : (a) normal ; (b) acute inflammation ; (c) subcute chronic 

inflammation ; (d) cirrhosis of the liver ; (e) nephrotic syndrome; (f) (J-myeloma; 

(g) Y-plasmocytoma[19]. 
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Figure 2.3.2 Stability and accuracy characteristics for implicit SUPG formulation. 
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Figure 2.3.3 Stability and accuracy characteristics for explicit 1-pass SUPG 
formulation. 
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Figure 2.3.4 Stability and accuracy characteristics for explicit 2-pass SUPG 
formulation. 
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13.6 Stability and accuracy characteristics 
formulation. 


for explicit 1-pass TW/SW 
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Figure 2.4.2 

Convection of a discontinuity . 












Figure 2.4.3 Initial condition for translating puff and rotating puff problems. 
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Figure 2.4.4 Translating puff elevation plots ( SUPG ). 
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Figure 2.4.5 Translating puff elevation plots ( TW/SW ). 
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Figure 2.4.6 Rotating puff elevation plots ( SUPG ). 
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Figure 2.4.7 Rotating puff elevation plots ( TW ). 
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Figure 2.4.8 Rotating puff elevation plots ( SW ). 
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WALL DETAILS 
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Figure 3.2. 1 Vorticity -stream function formulation. Wall boundary details. 



FINE MESH QOA- (ELEMENTS) 



Figure 3.5.1 Flow past a step. Problem description, boundary conditions, and mesh 
employed. 
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Re = 200 (coarse mesh) Steady state streamlines 



Corner streamlines 



Re — 200 (coarse mesh) Steady state iso — vorticity lines 



Figure 3.5.2 Flow past a step at Re = 200. Steady state streamlines and iso-vorticity lines 
(coarse mesh ). 
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Figure 3.5.3 Flow past a step at Re * 200. Streamlines at various time steps( fine mesh ). 
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Figure 3.5.3 Flow past a step at Re « 200. Streamlines at various time steps( fine mesh ). 
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Corner streamlines 



Figure 3.5.3 Row past a step at Re « 200. Streamlines at various time steps( fine mesh ). 
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Figure 3.5.3 Flow past a step at Re =* 200. Streamlines at various time steps( fine mesh ). 
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Re = 200 (fine mesh) Streamlines at t = 2.0 



Corner Streamlines 



- Re = 200 (fine mesh) Iso — vorticity lines at t = 2.0 



Figure 3.5.3 Flow past a step at Re « 200. Steady state streamlines and iso- vorticity lines. 
( fine mesh ). 
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DRIVEN CAVITY PROBLEM 



Figure 3.5.4 Driven cavity. Problem description and boundary conditions. 
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Figure 3.5.6 Driven cavity at Re * 400. Steady state streamlines and iso-vorticity lines. 




Figure 3.5.7 Driven cavity at Re = 800. Steady state streamlines and iso- voracity lines. 



Figure 3.5.8 Driven cavity at Re = 100 and 400. Steady state streamlines and iso-vorticity 
lines obtained by Hughes et aL [32]. 
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ZONE ELECTROPHORESIS 



Figure 4.1.1 Zone electrophoresis. A mixture S consisting of anions and cations to be 
separated is injected in the electrolyte consisting of P + and Q' ions. On the 
passage of an electric current the species separate out into zones depending on 
the charge and mobility. 
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MOVING BOUNDARY ELECTROPHORESIS 



Figure 4. 1 2 Moving boundary electrophoresis. The sample S consisting of species A and B 
is introduced into the cathode compartment. On the passage of current species A 
separates out fully because it has a higher mobility. This is followed by a 
mixture of A+B. 
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ISOTACHOPHORESIS 
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Figure 4.1.3 Isotachophoresis. The sample S consisting of A+B species is introduced 
between leading( L ) and tenninating( T ) anionic species. On the passage of 
current mixed zones are obtained( like in MBE). On achieving steady state 
species A and B separate out and move at the same constant speed. 
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ZONE ELECTROPHORESIS IN ONE DIMENSION 
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Figure 4.6.1 Zone electrophoresis in 1-D. Temporal evolution of species A* and B\ 
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ZONE ELECTROPHORESIS IN ONE DIMENSION 
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Figure 4.6.1 Zone electrophoresis in 1-D. Temporal evolution of P* and potential gradient 
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MOVING BOUNDARY ELECTROPHORESIS IN ONE DIMENSION 
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Figure 4.6.2 Moving boundary electrophoresis in 1-D. Temporal evolution of species A + 
and B*. 
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MOVING BOUNDARY ELECTROPHORESIS IN ONE DIMENSION 
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Figure 4.6 2 Moving boundary electrophoresis in 1-D. Temporal evolution of P 
and potential gradient 
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ISOTACHOPHORESIS IN ONE DIMENSION 



Figure 4.6.3 Isotachophoresis in 1-D. Initial condition showing concentration profiles of 
terminator T\ species B* and A'*’, and leader L + (in order of increasing 
mobility). 
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ISOTACHOPHORESIS IN ONE DIMENSION 
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Figure 4.6.3 Isotachophoresis in 1-D. Intermediate and final condition showing concentration 
profiles of terminator T\ species B + and A*, and leader L* (in order of 
increasing mobility). 
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Figure 4.6.4 Isotachophoresis in 1-D. Temporal evolution of potential gradient. 
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ISOTACHOPHORESIS IN ONE DIMENSION 




Figure 4.6.4 Isotachophoresis in 1-D. Temporal evolution of terminator T+ and leader L+. 
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ISOTACHOPHORESIS IN ONE DIMENSION 


Potential Gradient for ITP 



Figure 4.6.4 Isotachophoresis in 1-D. Temporal evolution of potential gradient 



OF POOR QUALITY 


ZONE ELECTROPHORESIS IN TWO DIMENSIONS 
Species A + Species B‘ 




Figure 4.6.5 Zone electrophoresis in 2-D. Temporal evolution of species A + and B . 
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Figure 4.6.5 Zone electrophoresis in 2-D. Temporal evolution of P* and potential gradient. 



Figure 4.6.6 Cross electrophoresis in 1*D. Temporal evolution of species A * and B\ 
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Figure 4.6.6 Cross electrophoresis in 1-D. Temporal evolution of P* and potential gradient 


